---
title: "Ideology and the Red Button"
author: "Onderco, M., Etienne, T., & Smetana, M."
date: "April 18, 2022"
output:
  html_document:
    data_print: paged
    toc: yes
    toc_depth: '3'
    toc_float: yes
  word_document: null
  html_notebook:
    number_sections: no
    theme: cerulean
    toc: yes
    toc_depth: 3
    toc_float: yes
  pdata_document:
    toc: yes
    toc_depth: '3'
    toc_float: yes
  latex_engine: xelatex
  pdf_document:
    toc: yes
    toc_depth: '3'
subtitle: Data and analyses for Foreign Policy Analysis
editor_options:
  markdown:
    wrap: 72
---

```{r, warning=F, echo=F, message=F}
# Script author: Tom W. Etienne

# Load packages
## reading
library(haven) # load spss files

## wrangling
library(dplyr) # data wrangling
library(tidyr) # data wrangling
library(tidyverse) # data wrangling

## analysis
library(psych) # cronbach's alpha

## output
library(ggplot2) # plots
library(sjPlot) # plot regression results
library(knitr) # tables


# turn off scientific notation
options(scipen=999)
options(digits=10)

# Set seed
set.seed(1)
```

\newpage

# Data loading and cleaning

```{r, warning=F}
data = as.data.frame(read_spss("./ONDERCO-ETIENNE-SMETANA_FPA_data_v6.sav"))
```

## Missing values

```{r, warning=F}
data$LR = ifelse(data$LR %in% 0:10, data$LR, NA)
```

## Add labels

```{r, warning=F}
# add a variable to the data with the labels of the LR-scale in 5 bins
LR_5bins_labels = as.data.frame(data$LR_5bins %>% attr('labels'))
data$LR_5bins_labels = row.names(LR_5bins_labels)[match(data$LR_5bins,LR_5bins_labels$`data$LR_5bins %>% attr("labels")`)]
```

\newpage

# Article

## Descriptives

```{r}
# sex
n_sex = as.data.frame(table(data$COUNTRY, data$sex_regr))
colnames(n_sex) = c("Country", "Sex", "n")
n_sex$Sex = ifelse(n_sex$Sex==0, "Male",
            ifelse(n_sex$Sex==1, "Female",NA))
kable(n_sex)

# age
(mean_age_DE = round(mean(data$age_cont[data$COUNTRY=="DE"]),0))
(sd_age_DE = round(sd(data$age_cont[data$COUNTRY=="DE"]),1))

(mean_age_NL = round(mean(data$age_cont[data$COUNTRY=="NL"]),0))
(sd_age_NL = round(sd(data$age_cont[data$COUNTRY=="NL"]),1))

# education
n_edu = as.data.frame(table(data$COUNTRY, data$edu_2_rec))
colnames(n_edu) = c("Country", "Education", "n")
n_edu$Education = ifelse(n_edu$Education==0, "No higher education",
                  ifelse(n_edu$Education==1, "Higher education",NA))
kable(n_edu)
```


## Cronbach's alpha

```{r}
items = c("scenario1", "scenario2", "scenario3", "scenario4")
psych::alpha(data[,names(data) %in% items])
```


## Table 2: DESCRIPTIVES FOR THE FOUR SCENARIOS OF NUCLEAR WEAPONS USE

The syntax for table 2 can be found in the code chunk below.

```{r, warning=F}
# scenario 1
table2.1 = data %>% filter(is.na(scenario1_bin)==F) %>%
  group_by(COUNTRY) %>% 
  count(scenario1_bin , wt = Combined_weight_trimmed_0.995) %>% 
  mutate(prop = n/sum(n)) %>%
  select(-n) %>%
  pivot_wider(names_from = c(COUNTRY, scenario1_bin), values_from = prop)

# scenario 2
table2.2 = data %>% filter(is.na(scenario2_bin)==F) %>%
  group_by(COUNTRY) %>% 
  count(scenario2_bin , wt = Combined_weight_trimmed_0.995) %>% 
  mutate(prop = n/sum(n)) %>%
  select(-n) %>%
  pivot_wider(names_from = c(COUNTRY, scenario2_bin), values_from = prop)

# scenario 3
table2.3 = data %>% filter(is.na(scenario3_bin)==F) %>%
  group_by(COUNTRY) %>% 
  count(scenario3_bin , wt = Combined_weight_trimmed_0.995) %>% 
  mutate(prop = n/sum(n)) %>%
  select(-n) %>%
  pivot_wider(names_from = c(COUNTRY, scenario3_bin), values_from = prop)

# scenario 4
table2.4 = data %>% filter(is.na(scenario4_bin)==F) %>%
  group_by(COUNTRY) %>% 
  count(scenario4_bin , wt = Combined_weight_trimmed_0.995) %>% 
  mutate(prop = n/sum(n)) %>%
  select(-n) %>%
  pivot_wider(names_from = c(COUNTRY, scenario4_bin), values_from = prop)

# combine
table2 = rbind(rbind(rbind(table2.1, table2.2), table2.3), table2.4)
table2$Scenarios = c("Scenario 1", "Scenario 2", "Scenario 3", "Scenario 4")
table2 = table2[,c(5,1,3,2,4)]
table2 = round(table2[,2:5],3)
kable(table2)
```


## Figure 1: WILLINGNESS TO USE NUCLEAR WEAPONS BY COUNTRY AND IDEOLOGY 

```{r, warning=F, fig.width=12}
# select data without missing values on the weight variable, and the plotted variable
data_selection = data[!is.na(data$LR_5bins_labels), ]

x_labels = c("Strongly \n disagree","Disagree","Tend to \n disagree", "Tend to \n agree", "Agree", "Strongly \n agree")

ggplot( data = data_selection, # data
        aes(x        = NW_willingness_index_avg, # x axis data
            col      = LR_5bins_labels, # visual grouping
            linetype = LR_5bins_labels, # different linetypes
            weights  = Combined_weight_trimmed_0.995)) + # weights
        geom_density(bw=.5, size=1) + # kernel density + line thickness
        labs(title = "Average willingness to use nuclear weapons",
             x = "",
             y = "Density", 
             col = "Left-right \nself-identification",
             linetype = "Left-right \nself-identification") + # labels legend title
        scale_color_manual   (values=c("#00A300", "#008b8b", "#000075", "#800000", "#ffa500")) + # line colors
        scale_linetype_manual(values=c(1,5,2,4,3)) + # linetypes
        guides(colour = guide_legend(override.aes = list(size=1, linetype=c(1,5,4,2,3)))) + # legend
        scale_x_discrete     (limits = x_labels) + # x axis labels
        facet_grid(cols=vars(COUNTRY)) + # two panels by country
        theme_bw()

```

## Table 3: RESULTS OF QUANTITATIVE ANALYSIS

The regressions are done in SPSS rather than in R. Please use the script ***FPA_regression_analyses.sps***.

## Figure 2: IDEOLOGY AND WILLINGNESS TO USE FORCE

```{r, warning=F, fig.width=12}
# create polynomial regressions to plot
poly1_DE = lm(data=data, weights = weight_trimmed_0.995_DE,
              NW_willingness_index_avg ~ age_cont+sex_regr+geo_DE_East
              +LR)
poly1_NL = lm(data=data, weights = weight_trimmed_0.995_NL,
              NW_willingness_index_avg ~ age_cont+sex_regr
              +LR)
poly3_DE = lm(data=data, weights = weight_trimmed_0.995_DE,
              NW_willingness_index_avg ~ age_cont+sex_regr+geo_DE_East
              +LR+eval(LR^2)+eval(LR^3))
poly3_NL = lm(data=data, weights = weight_trimmed_0.995_NL,
              NW_willingness_index_avg ~ age_cont+sex_regr
              +LR+eval(LR^2)+eval(LR^3))

# pull out predicted values
poly1_DE_preds = plot_model(poly1_DE, type = "pred", terms = c("LR"))
poly1_NL_preds = plot_model(poly1_NL, type = "pred", terms = c("LR"))
poly3_DE_preds = plot_model(poly3_DE, type = "pred", terms = c("LR"))
poly3_NL_preds = plot_model(poly3_NL, type = "pred", terms = c("LR"))

# calculate weighted means and add predicted values
avgs = data_selection %>% group_by(LR, COUNTRY) %>% 
  summarise(mean = weighted.mean(NW_willingness_index_avg, 
                                 w=Combined_weight_trimmed_0.995, na.rm=T))

avgs$poly1[avgs$COUNTRY=="DE"] = poly1_DE_preds$data$predicted
avgs$poly1[avgs$COUNTRY=="NL"] = poly1_NL_preds$data$predicted
avgs$poly3[avgs$COUNTRY=="DE"] = poly3_DE_preds$data$predicted
avgs$poly3[avgs$COUNTRY=="NL"] = poly3_NL_preds$data$predicted

# plot
ggplot(data=avgs, aes(x=LR)) +
  # data elements
  geom_point(aes(y=mean)) +
  geom_line (aes(y=poly1)) +
  geom_line (aes(y=poly3), lty=2) +
  facet_wrap(~COUNTRY) +
  # settings
  scale_x_continuous(breaks = seq(from=0,to=10,by=1)) +
  scale_y_continuous(breaks = seq(from=0,to=6,by=1),limits=c(0, 6)) +
  labs(title = "Willingness to use nuclear weapons",
       y = "Average willingness to use NW (index value)",
       x = "Left-right self-placement") +
  theme_bw() + theme(panel.grid.minor = element_blank())


```

\newpage

# Appendices

## Appendix Figure 1

```{r, warning=F, fig.width=12}
# scenario1
ggplot( data=data_selection, # data
        aes(x=scenario1, # x axis data
            #group=LR_5bins_labels, # grouping -> no longer needed?
            col=LR_5bins_labels, # visual grouping
            linetype=LR_5bins_labels, # different linetypes
            weights=Combined_weight_trimmed_0.995)) + # weights
  geom_density(bw=.5, size=1) + # kernel density + line thickness
  guides(colour = guide_legend(override.aes = list(size=1, linetype=c(1,5,4,2,3)))) + # legend line type & thickness
  labs(title = "Demonstrative scenario to de-escalate",
       x = "",
       y = "Density", 
       col = "Left-right \nself-identification",
       linetype = "Left-right \nself-identification") + # labels legend title
  scale_color_manual   (values=c("#00A300", "#008b8b", "#000075", "#800000", "#ffa500")) +
  scale_linetype_manual(values=c(1,5,2,4,3)) +
  scale_x_discrete     (limits=c("Strongly \n disagree","Disagree","Tend to \n disagree", "Tend to \n agree", "Agree", "Strongly \n agree")) + # x axis labels
  facet_grid(cols=vars(COUNTRY)) + # two panels by country
  theme_bw()
```

## Appendix Figure 2

```{r, warning=F, fig.width=12}
# scenario2
ggplot( data=data_selection, # data
        aes(x=scenario2, # x axis data
            #group=LR_5bins_labels, # grouping -> no longer needed?
            col=LR_5bins_labels, # visual grouping
            linetype=LR_5bins_labels, # different linetypes
            weights=Combined_weight_trimmed_0.995)) + # weights
  geom_density(bw=.5, size=1) + # kernel density + line thickness
  guides(colour = guide_legend(override.aes = list(size=1, linetype=c(1,5,4,2,3)))) + # legend line type & thickness
  labs(title = "Targeting Russian military units",
       x = "",
       y = "Density", 
       col = "Left-right \nself-identification",
       linetype = "Left-right \nself-identification") + # labels legend title
  scale_color_manual   (values=c("#00A300", "#008b8b", "#000075", "#800000", "#ffa500")) +
  scale_linetype_manual(values=c(1,5,2,4,3)) +
  scale_x_discrete     (limits=c("Strongly \n disagree","Disagree","Tend to \n disagree", "Tend to \n agree", "Agree", "Strongly \n agree")) + # x axis labels
  facet_grid(cols=vars(COUNTRY)) + # two panels by country
  theme_bw()
```

## Appendix Figure 3

```{r, warning=F, fig.width=12}
# scenario3
ggplot( data=data_selection, # data
        aes(x=scenario3, # x axis data
            #group=LR_5bins_labels, # grouping -> no longer needed?
            col=LR_5bins_labels, # visual grouping
            linetype=LR_5bins_labels, # different linetypes
            weights=Combined_weight_trimmed_0.995)) + # weights
  geom_density(bw=.5, size=1) + # kernel density + line thickness
  guides(colour = guide_legend(override.aes = list(size=1, linetype=c(1,5,4,2,3)))) + # legend line type & thickness
  labs(title = "Response to Russia's demonstration",
       x = "",
       y = "Density", 
       col = "Left-right \nself-identification",
       linetype = "Left-right \nself-identification") + # labels legend title
  scale_color_manual   (values=c("#00A300", "#008b8b", "#000075", "#800000", "#ffa500")) +
  scale_linetype_manual(values=c(1,5,2,4,3)) +
  scale_x_discrete     (limits=c("Strongly \n disagree","Disagree","Tend to \n disagree", "Tend to \n agree", "Agree", "Strongly \n agree")) + # x axis labels
  facet_grid(cols=vars(COUNTRY)) + # two panels by country
  theme_bw()
```

## Appendix Figure 4

```{r, warning=F, fig.width=12}
# scenario4
ggplot( data=data_selection, # data
        aes(x=scenario4, # x axis data
            #group=LR_5bins_labels, # grouping -> no longer needed?
            col=LR_5bins_labels, # visual grouping
            linetype=LR_5bins_labels, # different linetypes
            weights=Combined_weight_trimmed_0.995)) + # weights
  geom_density(bw=.5, size=1) + # kernel density + line thickness
  guides(colour = guide_legend(override.aes = list(size=1, linetype=c(1,5,4,2,3)))) + # legend line type & thickness
  labs(title = "Targeting Kaliningrad after Russian NW use",
       x = "",
       y = "Density", 
       col = "Left-right \nself-identification",
       linetype = "Left-right \nself-identification") + # labels legend title
  scale_color_manual   (values=c("#00A300", "#008b8b", "#000075", "#800000", "#ffa500")) +
  scale_linetype_manual(values=c(1,5,2,4,3)) +
  scale_x_discrete     (limits=c("Strongly \n disagree","Disagree","Tend to \n disagree", "Tend to \n agree", "Agree", "Strongly \n agree")) + # x axis labels
  facet_grid(cols=vars(COUNTRY)) + # two panels by country
  theme_bw()
```

## Appendix Figure 5

```{r, warning=F}

# Label vote intent
data$DE_voteintent =  ifelse(data$vote_intent_DE == 1, "CDU", 
                      ifelse(data$vote_intent_DE == 2, "SPD", 
                      ifelse(data$vote_intent_DE == 3, "AfD", 
                      ifelse(data$vote_intent_DE == 4, "Linke", 
                      ifelse(data$vote_intent_DE == 5, "FDP", 
                      ifelse(data$vote_intent_DE == 6, "Die Grunen", 
                      ifelse(data$vote_intent_DE == 11, "CSU",NA)))))))

# select data
data_selection = data[!is.na(data$DE_voteintent),]

## make plot
ggplot(data_selection, aes(x=DE_voteintent, 
                  y=LR)) + # data and variables # weight=Combined_weight_trimmed_0.995
  ggtitle("Left-right self placement by vote intention (unweighted)") +
  scale_y_continuous(name="Left-right self-placement", breaks=seq(0,10,1),limit=c(0,10)) + # y axis
  scale_x_discrete(name ="Vote intention", 
                    limits=c("Linke","SPD","Die Grunen","CDU","FDP","CSU","AfD")) + # x axis
  geom_violin(scale="area", bw=.4) + # violing with smoothing kernel bandwidth, and area equal
  geom_boxplot(width = .1, outlier.shape = NA) + # boxplot no outliers shown
  theme_bw()


```

## Appendix Figure 6

```{r, warning=F, warning=F}

data$NL_voteintent =  ifelse(data$vote_intent_NL == 1, "VVD",
                      ifelse(data$vote_intent_NL == 5, "GL",
                      ifelse(data$vote_intent_NL == 8, "CU",
                      ifelse(data$vote_intent_NL == 3, "CDA",
                      ifelse(data$vote_intent_NL == 6, "SP",
                      ifelse(data$vote_intent_NL == 7, "PvdA",
                      ifelse(data$vote_intent_NL == 2, "PVV",
                      ifelse(data$vote_intent_NL == 13, "FvD",
                      ifelse(data$vote_intent_NL == 4, "D66",
                      ifelse(data$vote_intent_NL == 9, "PvdD",
                      ifelse(data$vote_intent_NL == 11, "SGP",
                      ifelse(data$vote_intent_NL == 10, "50PLUS",NA))))))))))))

# select data
data_selection = data[!is.na(data$NL_voteintent),]

## make plot
ggplot(data_selection, aes(x=NL_voteintent, 
                  y=LR)) + # data and variables # weight=Combined_weight_trimmed_0.995
  ggtitle("Left-right self placement by vote intention (unweighted)") +
  scale_y_continuous(name="Left-right self-placement", breaks=seq(0,10,1),limit=c(0,10)) + # y axis
  scale_x_discrete(name ="Vote intention", 
                    limits=c("GL","SP","PvdD","PvdA","D66","CU","50PLUS","CDA","VVD","FvD","SGP","PVV")) + # x axis
  geom_violin(scale="area", bw=.4) + # violing with smoothing kernel bandwidth, and area equal
  geom_boxplot(width = .1, outlier.shape = NA) + # boxplot no outliers shown
  theme_bw()
```

## Appendix 1

The code for the regressions can be found in ***FPA_appendix_regression.sps***.

## Appendix 2

The code for the regressions can be found in ***FPA_appendix_regression.sps***.

## Appendix 3

The code for the regressions can be found in ***FPA_appendix_regression.sps***.